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Abstract 

Implementing  large  scale  LP  algorithms  is  a difficult,  laborious 
and  often  poorly  rewarded  task.  This  is  particularly  true  of  algorithms 
which  exploit  the  structure  of  the  LP  matrix.  For  this  reason  many 
algorithms  have  been  proposed  in  the  literature,  but  few  have  been 
turned  into  good  computer  codes.  Very  little  is  known  about  the  rela- 
tive performance  of  different  algorithms.  In  this  paper  we  discuss 
some  of  the  suggestions  that  have  been  made  for  alleviating  this  problem 
and  describe  an  approach  based  upon  a carefully  defined  collection  of 
subroutines  which  are  designed  to  aid  the  task  of  implementing  and 
comparing  LP  algorithms.  These  subroutines  or  modules  may  be  regarded 
as  the  'primitives’  of  a language  for  Implementing  experimental  LP 
algorithms,  particularly  algorithms  which  exploit  matrix  structure. 

A set  of  such  modules  is  described.  These  have  been  implemented  in 
FORTRAN,  and  user  documentation  is  available. 
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MODULES  TO  AID  THE  IMPLEMENTATION  OF  LP  ALGORITHMS 


by 

L.  Nazareth 


1.  Introduction 

Implementing  large  scale  LP  algorithms  Is  a difficult,  laborious 
and  often  poorly  rewarded  task.  This  is  particularly  true  of  algorithms  which 
exploit  the  structure  of  the  LP  matrix.  For  this  reason  many  algorithms  have 
been  proposed  in  the  literature,  but  few  have  been  turned  into  good  computer 
codes.  Very  little  is  known  about  the  relative  performance  of  different 
algorithms. 

In  this  paper  we  discuss  some  of  the  suggestions  that  have  been  made 
for  alleviating  this  problem,  and  describe  an  approach  based  upon  a carefully 
defined  collection  of  subroutines  which  are  designed  to  aid  the  task  of 
implementing  and  comparing  LP  algorithms.  These  subroutines  or  modules  may 
be  regarded  as  the  'primitives'  of  a language  for  implementing  experimental 
LP  algorithms,  particularly  algorithms  which  exploit  matrix  structure.  A set 
of  such  modules  is  described.  These  have  been  implemented  in  FORTRAN,  and 
user  documentation  is  provided.  (We  developed  these  modules  as  a first  stage 
in  the  implementation  of  a nested  decomposition  code  for  staircase  structures, 
but  they  were  purposely  designed  for  use  in  other  contexts,  and  they  may  also 
be  a useful  educational  aid.) 

This  has  been  primarily  a software  organization  and  development  effort 
and  we  have  drawn  upon  the  work  of  many  different  researchers  in  the  field. 

We  do  not  claim  that  the  modules  are  quality  software,  but.  we  hope  that  they 
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will  evolve  in  the  direction  of  quality  software,  e.g.,  through  much  more 
exhaustive  testing  of  the  modules. 

In  the  concluding  section,  we  discuss  future  directions  and  the 
benefits  to  be  gained  from  an  effort  such  as  this  one. 


2 . Background 

We  have  discussed  the  overall  framework  for  optimization  software 
development  elsewhere,  see  Nazareth  [1],  [2],  In  summary,  we  distinguish 
between  implementations  designed  primarily  to  aid  algorithms  or  code  develop-  , 

ment  and  implementations  which  are  intended  for  production  runs  and  are  thus 
designed  primarily  to  solve  user  problems.  The  former  are  called  algorithm/ code 
oriented  implementations  (or  sometimes  experimental  implementations)  and  the 
later  user/problem  oriented  implementations.  Each  of  these  can  be  further 
tailored  to  a particular  compiler  and  machine  configuration,  resulting  in 
different  program  realizations  (see  Boyle  and  Dritz  [3]). 

In  Nazareth  fl } , three  different  approaches  to  algorithm/code  oriented 
software  were  also  discussed.  In  the  first  approach  one  seeks  to  develop  a 
suitable  high  level  language.  This  language  permits  highly  readable  programs 
to  be  written  with  relative  ease  in  the  vernacular  of  applied  mathematics, 
and  serves  as  a medium  for  communicating  algorithmic  ideas  precisely  (MPL, 
see  Dantzig  et  al.  [4],  is  an  example  of  such  a language.  Another  example 
is  the  Speakeasy  language  of  Cohen  et  al.  [5],  This  latter  language  also 
provides  a mechanism  which  enables  a user  to  extend  the  language  to  suit 
his  individual  needs).  In  the  early  creative  stages  of  algorithm  development,  j 

such  a language  is  very  useful,  since  new  ideas  can  be  quickly  implemented 
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and  tested  out, and  the  computational  experience  obtained  often  results  in  new 
Insights  and  developments.  Thus  the  major  features  of  an  algorithm  can  be 
laid  out.  Note  however  that  a 'quick  and  dirty'  implementation  can  usually 
be  run  only  on  toy  problems,  and  that  procedures  which  are  numerically  sound 
and  efficient  in  terms  of  time  and  storage,  are  difficult  to  write  in  any 
language.  Other  pros  and  cons  of  this  first  approach  are  discussed  in 
Nazareth  [1]. 


In  the  second  approach,  one  attempts  to  introduce  some  limited 
flexibility  into  software  designed  primarily  to  solve  user  problems,  by  making 
some  of  the  high  level  components  available  to  the  algorithm  developer.  For 
example,  in  the  MPSX/370  system,  a number  of  high  level  routines,  listed  in 
Figure  1,  can  be  accessed  through  a control  language  and  used  to  construct 
algorithms.  Another  example  is  the  SEXOP  system,  Madsen  [6],  see  Fig.  2,  in 
which  the  major  components  are  geared  toward  problem  manipulation.  Note  that 
in  both  cases  one  is  greatly  constrained  by  the  overall  system,  particularly 
by  its  data  structures. 


The  third  approach  complements  the  first,  by  placing  much  more 
emphasis  on  aids  for  developing  a numerically  sound  experimental  implemen- 
tation, once  some  of  the  major  features  of  an  algorithm  have  been  laid  out  (for 
example,  using  an  implementation  in  MPL) . In  this  third  approach,  one 
seeks  to  Identify  the  components  that  are  used  to  build  an  LP  algorithm, 
to  specify  them  cleanly  and  carefully,  and  to  implement  them  in  a manner 
which  makes  them  flexible  and  easy  to  use.  These  modules  can  also  be 
thought  of  as  the  'primitives'  or  'basic  operators'  of  a language  for  build- 
ing LP  algorithms.  Since  our  aim  is  to  produce  software  which  is  immediately 
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FIGURE  2:  MAJOR  COMPONENTS  OF  SEXOP 


^SETUP 

^ADDCOL 

✓REPCOL 

SETCOL 

XCOST 

SCRHS 

'LINEAR 

PERTURB 

TRADOF 


(read  in  LP) 

(add  new  colunm  ) 
(replace  colunm  ) 
(examine  column  ) 


(change  objective  function  ) 
(change  rhs) — " 


(Involve  solver) 

(parametric  analysis  of  rhs) 


(parametric  analysis  of  objective) 


and  widely  usable,  we  think  that  they  should  be  implemented  in  FORTRAN 
(they  could  equally  well  be  implemented  in  another  high  level  language  and 
our  effort  can  be  regarded  as  another  contribution  to  identifying  what 
should  go  into  a language  for  optimization;  see-  also  Land  and  Powell  [7], 
Dantzig  et  al.  [8]).  Experimental  codes  which  are  constructed  from  these 
modules  should  be  executable  on  real  life  problems,  not  just  on  toy  data; 
thus  attention  should  be  paid  both  to  efficient  data  representations  and  to 
robustness.  Finally  the  modules  should  serve  as  an  educational  aid. 

A set  of  such  modules  are  described  in  the  remainder  of  this  paper. 

How  the  set  can  be  enlarged,  and  ways  to  make  them  easier  to  use,  are  discussed 
in  Section  5. 


3.  Description  of  Modules 

The  modules  fall  into  three  categories  as  described  below.  The 
following  overview  is  brief  and  assumes  that  the  reader  is  familiar  with 
the  theory  of  LP.  For  more  detail  on  each  module  consult  the  documentation, 
and  for  a surfeit  of  detail  check  the  listings  (see  Section  4 ) . 

3.1.  Problem  Oriented  Modules 

Problem  oriented  modules  serve  as  an  'interface'  between  the  user's 
LP  problem  and  the  routines  which  solve  it. 

LP  problems  are  usually  specified  on  a tape  in  standard  MPS  input 
format  (see  [9]  and  Figure  4)  while  LP  routines  usually  work  on  packed  matrix 
representations.  A set  of  modules  have  been  provided  to  aid  in  the  reading 
of  an  LP  tape  and  setting  up  of  the  data  structures. 
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Typical  structured  LP  matrices  are  illustrated  in  Fig.  3. 
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FIGURE  3 


An  LP  routine  designed  to  take  advantage  of  structure  may  require 
a particular  representation  of  data.  Therefore  it  would  be  inappropriate 
to  provide  a general  input  routine.  Instead  we  provide  components  from 
which  a suitable  input  routine  can  be  built.  These  modules  can  be  used  in 
conjunction  with  those  of  Section  3.2.1,  to  set  up  the  data  representation 
as  required.  For  example,  given  a block  angular  LP  matrix,  it  may  be  necessary 


to  set  up  each 


C) 


as  an  LP  matrix  with  consecutively  numbered  rows.  A call  to  module  PREADC  (see 
below)  will  set  up  a packed  matrix  with  rows  numbered  according  to  the  original 
matrix.  Following  this  by  a call  to  ADRNDX  (see  Section  3.2.1)  will  renumber 
the  rows  to  be  consecutive. 


The  problem  oriented  modules  In  this  version  of  the  collection  are 
as  follows:  (They  were  adapted  In  major  part  from  the  input  routines  of 
MINOS,  Saunders  [10].  However,  since  we  segmented  the  routine  and  altered 
several  portions  of  the  code  to  suit  our  needs,  responsibility  for  errors 
rests  with  us,  and  should  in  no  way  reflect  upon  the  source  of  the  code.) 


o PREADR 
o PREADC 

o PRDRHS 
o PREADB 

o PCHKST 


Read  the  ROWS  section  of  an  LP  matrix 

Given  a specified  set  of  columns  in  the  COLUMNS  section,  build 
a column  llst/row  index  packed  data  structure. 

Read  the  RHS  section  of  an  LP  matrix 

Read  in  the  BOUNDS  section  of  an  LP  matrix  (see  also  Section 
3.2.2). 

Perform  checks  and  gather  statistics  about  the  matrix. 


The  manner  in  which  a matrix  is  packed  is  illustrated  in  Fig.  3 for 
a specific  example.  By  writing  an  input  routine  which  calls  the  problem 
oriented  modules  the  user  can  translate  the  MPS  input  into  a packed  data 
representation,  suitable  for  his  LP  algorithm. 


3.2.  Algorithm  Oriented  Modules 

These  are  the  basic  building  blocks  from  which  algorithms  are 
constructed.  They  fall  into  two  categories. 
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FIGURE  4 


min  X1  + x2  + x3 


■^^Column  Names 
Row  name  8^''-'-^ 


TABLEAU 


CLM1  CLM2  CLM3  RTH 


s.t.  2x,  + 3x_  < 10  i 

1 3 _ ; 

OBJ  . 

1. 

1. 

1. 

0 

4x2  + 5x3  < 20  ! 

RWN1 

2. 

0 

3. 

< 

10. 

100  > x,  > 0,  x.  > 0 ’ 

— 1 — 2 — i 

RWN2 

0 

4. 

5. 

_< 

O 

CN 

Example  of  Sample  MPS  Input 


(further  details  are  given  in 
Chapter  I of  the  documentation) 

NAME  LP 
ROWS 


L RWN1 

a 

L RVKW  ^ 

mi  TTMMQ 

a+2 

a+4 

CLM1 

OBJ 

1.0 

RWN1 

2.0 

CLM2 

OBJ 

1.0 

RWN2 

4.0 

a+7 

CLM3 

OBJ 

1.0 

RWN1 

3.0 

HE 

CLM3 

RWN2 

5.0 

RTH  RWN1  10.0 

RTH  RWN2  20.0 

BOUNDS 

UP  B”N  CLM1  luO. 

ENDATA 


Packed  Representation  of  Above  Matrix, 

Excluding  RHS  (column  list/row  index  data 
structure)  (further  details  are  given  in 
Chapter  I of  the  documentation) 

Column  Pointers  Matrix  Elements  Row  Indices 
I 1 a >i 1 a -*-i 1 


1. 

a -► 

i 

2. 

2 

1. 

1 

4. 

3 

1. 

1 

3. 

2 

5. 

3 

~~ *A 

HA 

Pointer  to  first 
unused  element  \ 
of  this  array 


Thus  the  third  column  (called  CLM3)  starts 
at  element  a+4  of  the  array  called  'matrix 
elements.'  This  column  has  three  elements 
whose  corresponding  indices  are  given  by  the 
elements  of  the  array  called  'row  indices.' 
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3.2.1.  Data  Structure  Manipulation  Modules 
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These  carry  out  a number  of  operations  on  packed  LP  matrices.  By 
using  them  It  is  possible  to  make  a distinction  between  devising  a strategy 
and  the  actual  mechanics  of  implementing  it. 

The  routines  that  are  included  in  the  current  version  of  the  collection 
are  as  follows: 


ADCONC  — 


ADRNDX  — 
ADTNTF  — 


ADUPKC  — 
ADDELC  — 


concatenate  two  packed  data  structures  and  return  result 
in  first  one 

reindex  rows  in  a packed  data  structure 

convert  packed  data  structure  to  element/row  index/ column 
index  data  structure 

unpack  a column  of  a packed  data  structure 
delete  a coluun  in  a packed  data  structure. 


Obviously  there  are  many  other  operations  of  this  sort.  The  above 
are  some  which  arise  in  the  implementation  of  decomposition  algorithms  (c.f. 
Introduction) . 


3.2.2.  Simplex  Modules 

Algorithms  for  structured  LP  are  usually  based  upon  repeated 
calls  to  a routine  which  performs  one  or  more  Iterations  of  the  revised 
simplex  method.  This  routine  must  be  specifically  designed  to  meet  the  needs 
of  the  calling  algorithm  for  structured  LP,  but  the  task  of  implementing  it 
is  made  natch  easier  if  the  following  set  of  basic  operations  are  available: 
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• MODRHS 

• FORMC 

• PRICE 

• CHVZR 

• UPBETA 


given  values  of  non-basics,  develop  modified  right-hand  side, 
form  cost  vector  for  Phases  1 or  2 of  Simplex  method 
price  out  columns 

determine  which  column  leaves  basis 
update  current  approximation  to  solution 


All  these  modules  utilize  a common  data  structure,  as  shown  in  Fig.  6,  and 
it  is  assumed  that  the  packed  LP  matrix  is  in  'computational  canonical  form,' 
as  explained  in  Fig.  6.  In  particular: 


1.  All  bounds  on  the  structural  variables  x are  set  in  BL  and  BU.  If  a 
variable  is  unbounded  above  or  below,  the  corresponding  element  of  BL  or 
BU  is  set  to  a machine  representation  of  + «•.  (This  is  done  by  PREADB.) 

2.  A full  identity  matrix  for  the  logical  variables  (z^,^)  is  assumed  to  be 
written  at  the  start  of  the  matrix.  The  bounds  on  these  variables  are 
determined  by  the  type  of  rows  (as  specified  in  the  rows  section  and 

read  by  PREADR) . Again  see  Fig.  5.  Thus  no  distinction  is  really  made 
between  positive  and  negative  slacks  and  artificial  variables.  They  simply 
have  different  bounds  that  they  must  satisfy. 

An  extension  of  the  simplex  method  is  inherent  in  the  design  of  the 
modules  and  the  associated  data  structure,  in  that  non-baslc  variables  can  be 
fixed  at  values  between  their  bounds,  as  specified  by  PEG.  This  is  related 
to  the  superbasic  variables  of  Murtagh  and  Saunders  [11],  but  the  latter  are 
used  in  a more  powerful  way,  since  an  optimization  is  carried  out  in  the  sub* 
space  they  define.  The  use  of  pegged  variables  Involves  some  straightforward 
extensions  to  PRICE,  CHUZR  and  UPBETA  (see  documentation  for  more  detail) . 

PEG  contains  the  current  value  of  every  variable  in  the  problem,  both  loglcals 
and  structurals.  Thus  there  is  some  redundancy  of  information  stored;  but 
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FIGURE  5:  COMPUTATIONAL  CANONICAL  FORM 
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FIGURE  ft.  DATA  STRUCTURE  FOR  SIMPLEX  MODULES 


this  is  not  too  great  a penalty  to  pay  in  our  experimental  system,  given  the 
added  flexibility  that  PEG  makes  possible. 


2.2.2.  Sparse  Linear  Algebra 

The  routines  of  Reid  [12],  implementing  LU  factorization  of  sparse 
matrices  and  Bartel s-Golub  updating  are  very  suitable  for  our  purposes,  and  we 
only  provide  an  interface  for  converting  a basis  in  our  representation  to  one 
needed  by  Reid's  routine  LA05AD. 


A.  Documentation  and  Listings 

In  addition  to  this  paper,  machine-readable  documentation  has  been 
written,  organized  as  follows: 


Chapter  1:  User  documentation  giving  briefly  the  purpose,  usage  and  algorithmic 
details  for  each  module.  Each  group  of  modules  is  preceded  by  an 
Introductory  section  giving  background  material. 

Chapter  2:  Coding  and  documentation  conventions. 


Chapter  3:  Testing  programs  and  output  they  produce. 


Chapter  A:  Listings  of  modules. 
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Future  Directions  and  Conclusions 


One  obvious  direction  for  further  development  Is  the  addition  of  new 
v modules,  for  example,  problem  oriented  modules  to  output  solutions,  more 

extensive  data  manipulation  modules  (e.g.  to  reorder  a matrix  by  permutations), 
additional  simplex  modules  (e.g.  to  deal  with  dual  methods),  additional  sparse 
linear  algebra  routines,  and  so  on. 

A second  direction  for  expansion  is  to  make  the  modules  easier  to  use, 
by  deploying  them  within  a suitable  host  system  which  takes  over  much  of  the 
bookkeeping  tasks  of  a module  and  shields  the  user  from  them,  e.g.  the  Speakeasy 
system  of  Cohen  et  al.  [5].  In  addition,  the  host  system  provides  many  additional 
facilities,  e.g.  graphical  aids,  extensive  matrix  manipulations,  etc.  In  the 
calling  sequence  of  a typical  FORTRAN  module,  many  parameters  correspond  to 
work  vectors,  work  arrays,  dimensioning  information,  switches  and  error  flags. 

The  Speakeasy  language  provides  a very  convenient  mechanism  for  writing  an 
interface  to  a FORTRAN  subroutine.  In  essence  the  call  can  be  redefined  and 
often  greatly  simplified,  since  tasks  of  allocating  work  storage,  parameter 
checking,  etc.,  can  be  assumed  by  the  interface.  The  MPL  language  [4]  is 
another  example  of  a suitable  host  environment. 

The  effort  described  in  this  paper  is  limited  in  scope  and  suffers  from 
many  shortcomings.  However  we  would  like,  in  conclusion,  to  emphasize  three 
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benefits  that  are  gained  from  mathematical  software  patterned  along  these  lines: 
a)  It  is  of  use  to  investigators  in  Systems  Optimization  Laboratories  who 
wish  to  implement  optimization  algorithms  and  study  their  behavior.  In 
particular  we  plan  to  use  our  modules  to  implement  algorithms  for  structured 
LP  based  upon  the  Dantzlg-Wolfe  Decomposition  principle. 


b)  It  is  a useful  educational  aid;  in  particular  our  modules  serve  as  a mini- 

f. 

tutorial  on  the  implementation  of  LP  algorithms,  and  we  have  already  used 


them  for  this  purpose  in  the  lecture  hall, 
c)  The  discipline  of  systematizing  and  preparing  codes  raises  many  questions, 


which  in  turn  spurs  further  research.  The  PEG  array,  see  Fig.  6,  is  a small 


example  of  this.  This  point  is  also  discussed,  within  the  context  of  the 


UNPACK  project,  by  Stewart  [13,  p.  7]. 
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